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We examine the quantum entanglement spectra and Wannier functions of the square lattice Hofs- 
tadter model. Consistent with previous work on entanglement spectra of topological band structures, 
we find that the entanglement levels exhibit a spectral flow similar to that of the full system's energy 
spectrum. While the energy spectra are continuous, with open boundary conditions the entangle- 
ment spectra exhibit discontinuities associated with the passage of an energy edge state through 
the Fermi level. We show how the entanglement spectrum can be understood by examining the 
band projectors of the full system and their behavior under adiabatic pumping. In so doing we 
make connections with the original TKNN worP on topological two-dimensional band structures 
and their Chern numbers. Finally we consider Wannier states and their adiabatic flows, and draw 
connections to the entanglement properties. 



1. INTRODUCTION 

In the presence of a uniform magnetic field, the en- 
ergy spectrum of a noninteracting two-dimensional elec- 
tron gas is arranged in discrete, equally spaced Lan- 
dau levels. The Hall conductivity of n filled Landau 
levels is a xy = ne 2 /h. A discretized version of this 
model, due to Hofstadter^, has electrons hopping on a 
two-dimensional lattice with complex tight-binding am- 
plitudes tij — t e lAij , such that the magnetic flux through 
each unit cell is a rational multiple p/q of the Dirac 
flux quantum. The gauge field Aij can be made peri- 
odic by choosing a magnetic unit cell comprising q struc- 
tural unit cells of the lattice. For a lattice with an r 
element basis, this results in qr energy subbands which 
in general do not cross, a consequence of the Wigner- 
von Neumann theorem; the continuum limit is recovered 
at low energies for q — > oo. Plotting these energies as 
a function of <p = 2np/q yields the famous 'Hofstadter 
butterfly'. As shown in a seminal paper by Thouless et 
al^ (TKNN), to each band index j there corresponds 
an integer Chern number Cj , which physically represents 
the contribution to the Hall conductivity when band j 
is filled. The main differences with respect to the con- 
tinuum arc (i) the tight-binding subbands are dispersive, 
and (ii) whereas Cj = 1 for each Landau level in the 
continuum, the Chern indices of the TKNN bands are in 
general nonuniform. 

The Chern number is an integer invariant which re- 
flects aspects of the bulk band topology. As such it is 
robust and invariant with respect to parameter varia- 
tions which do not collapse the band gaps. The nontrivial 
bulk topology is also manifested at the edge. Hatsugai 3 
showed that the number of edge modes interpolating be- 
tween bulk bands separated by a gap is equal to the sum 
of the Chern indices of all bands below that gap. The 
spectral flow of the edge energy levels as a function of 
the momentum parallel to the edge is also reflected in 
the behavior of the quantum entanglement spectrurrPEl 
of the many-body reduced density matrix obtained by 
partitioning the system along a translationally-invariant 



boundary. For noninteracting fermions, the spectrum 
of the reduced density matrix itself corresponds to that 
of a noninteracting 'entanglement Hamiltonian' deter- 
mined by the one-body correlation matrix of the origi- 
nal systerrPED However there are also exceptions to the 
edge-entanglement correspondence. For example, the en- 
tanglement spectrum has protected midgap modes for a 
system wi th inv ersion symmetry even if the edge modes 
are gappecP^M In certain cases, one also has to tune the 
boundary conditions for a system with nontrivial topol- 
ogy in order for its energy edge modes to be gapless^, 
while such tuning is not required to observe the entangle- 
ment spectral flow. Thus in certain sense, the entangle- 
ment spectrum is a more robust test of the bulk topology. 

Once one specifies the wavevector k± along the trans- 
lationally invariant partition boundary, the entanglement 
Hamiltonian becomes effectively onc-dimcnsional, and 
the localization properties of such states can be con- 
sidered from the perspective of a Wannier basis^, and 
several recent studie s of t opological insulators have in- 
voked Wannier states^^ES! in their analyses. While non- 
vanishing Chern numbers provide an obstruction which 
rules out exponentially localized Wannier states in higher 
dimensions-^, the one-dimensional entanglement eigen- 
states at fixed k± can be so localized, and for topologi- 
cally nontrivial bulk band topologies, their Wannier cen- 
ters exhibit a spectral flow similar to that observed in 
the edge and entanglement spectra. 

In this paper, we investigate the spectral flow of entan- 
glement levels and Wannier states derived from the two- 
dimensional square lattice Hofstadter model. We identify 
the correspondence between these flows and the Chern 
numbers of the bulk bands, and investigate the effect 
of adiabatic pumping on the wavefunctions of different 
energy bands of the system. We consider both energy 
eigenstates as well as 'entanglement eigenstates' of the 
corresponding reduced density matrix which results from 
a spatial partitioning of the system into two parts. 
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2. HOFSTADTER MODEL AND ITS 
ENTANGLEMENT SPECTRUM 

2.1. Hofstadter model 

The Hofstadter modeP is a discrete model of electrons 
in two space dimensions and in the presence of magnetic 
flux. It is defined by a lattice tight-binding Hamiltonian, 

(«•) 

where = |tjj| is the complex hopping amplitude 
between sites i and j. The U(l) flux <\> p through a pla- 
quette p is the product Yl dp e lAij over a counterclockwise 
path of links along its boundary, dp. We shall only con- 
sider the case of uniform amplitude hopping, i.e. = I 
In the continuum limit, energy eigenstates of ballis- 
tic electrons collapse into macroscopically degenerate, 
equally spaced Landau levels. The degeneracy of each 
Landau level is N L = BQ/(/) Q , where B is the mag- 
netic field, f2 the total area covered by the system, and 
4> Q = hc/e is the Dirac flux quantum. The spectral flow 
of entanglement eigenstates in this limit was investigated 
by Rodriguez and Sierra 1 -. On the lattice, the Landau 
levels are no longer degenerate, but form magnetic sub- 



bands, each subband accommodating N L states. The 
model may be defined on any lattice, but for definiteness 
we consider the square lattice. Our principal results do 
not depend qualitatively on the underlying lattice struc- 
ture. (See Appendix [B] for the case of triangular lattice.) 

As is well-known, while 4> p is periodic on the scale of 
the structural unit cell, the vector potential A^ is not. 
However, if the flux <j> per plaquette is uniform and is 
2ir times a rational number p/q, a gauge can be chosen 
where A^j is periodic on the scale of a 'magnetic unit cell' 
comprising q elementary structural cells. For example, 
one can choose 

Aj = <f>yiS XitXj+ i 5 VitVj , (2) 

where {Xi,Vi) are integer coordinates for lattice site i. 
The magnetic unit cell is then a 1 x q tower of lattice cells, 
and one obtains a q x q Hamiltonian matrix with nonzero 
matrix elements H nn = — 2cos(fc a; + n<p), H Uyn+ i = — I, 
iiiv,i = —e lky , and remaining elements determined by 
hermiticity. This results in q magnetic subbands with 
dispersion e a {k x , k y ). Here we are concerned with entan- 
glement spectra, and to this end we consider a cylinder 
with periodic boundary conditions in the x-direction and 
N y sites in the y-direction. The Hamiltonian matrix is 



H(k x ,N y ,z) 



( 2 cos(k x 
I 





I 

2cos(fc a; + 2(j>) 1 

1 





(3) 



I 2 cos(fc x + N y (f>)J 



Here, k x is the Bloch phase in the x direction, N y is the 
number of lattice sites in the y direction, and z controls 
the boundary condition on y: z = for open boundary 
conditions, z = 1 with N y mod q — for periodic bound- 
ary conditions; unimodular complex z can be interpreted 
as flux threading the compactified cylinder. Note that 
for z € K the Hamiltonian is real and hence the eigen- 
functions may be chosen to be real as well. 

In what follows, we focus on the case with integer 
number K y of magnetic unit cells in the y direction 
(N y = qK y ), which is more convenient for switching be- 
tween open and periodic boundary. In ref. |3j the edge 
solution is derived with the requirement of commensu- 
rability, i.e., N y — qK y — 1, in order to exploit some 
structure in the transfer matrix formalism. This restric- 
tion can be lifted in the thermodynamic limit N y — > oo, 
where the spectrum of edge states localized at y = 1 is 
unchanged, while that of states localized at y = N y is 



shifted in k x . A proof is given in Appendix [A] 



2.2. Entanglement spectrum 

Imagine partitioning the sites of our lattice into two 
groups, A and B. If a many-body wavefunction |^) can 
be written as a direct product | V I' A ) ® \^ B ), the wave- 
function is said to be unentangled with respect to this 
partition. More generally, let g = | V I')(4'| be the pro- 
jector onto |$). Tracing out over the B sites yields the 
reduced density matrix g = Tr B g, whose eigenvalues con- 
stitute the entanglement spectrunM The von Neumann 
entropy, S = — Trplng, provides a measure of the degree 
of entanglement. 

If the boundary between A and B is irregular, then 
translational invariance is completely broken, but if the 
boundary is such that the A region remains periodic in 
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FIG. 1: (Color online) Energy levels (top) and entanglement 
occupancies (bottom) for the square lattice Hofstadter model 
with flux p/q = 3/7 per plaquette, on a cylinder of height 
N y = 28 as a function of the conserved crystal momentum 
k x . The Fermi level Ef lies inside the third gap (grey line), 
and the occupied bands below Ef contribute a total Chern 
number C occ = 1. The energy spectrum (top) is shown for 
open boundary conditions, with the black dots indicating bulk 
levels, the red lines indicating edge levels localized along the 
lower edge (y — 1), and the blue lines indicating edge levels 
localized along the upper edge (y = N y — 28). The vertical 
gray line marks the value of k x where the lower edge mode 
crosses the Fermi level. The entanglement occupancies f a 
(bottom) are computed for the lower half of the cylinder (1 < 
V < 14), color and symbol-coded according to a. Although 
the overall flow appears continuous, there is a discontinuity 
in the occupancies f a (k x ) where the lower edge mode crosses 
Ef, resulting in a sudden color change in the plot. 



one direction (or more than one, in the case of systems in 
more than two space dimensions) , then the eigenstates of 
g can be classified by a corresponding crystal momentum, 
and one can investigate the spectral flow of the entangle- 
ment levels^. 



FIG. 2: (Color online) Same as in Fig. [T] but with N y = 29 
and M = 15. 



2.2.1. The correlation matrix method 



A general recipe for computing the reduced density 
matrix for noninteracting Fermi systems has been de- 
rived by Cheong and Henley^. Let / and J denote sites 
in the full system, whose Hamiltonian is H — Hjj^Cj. 
The statistics of this Hamiltonian are then completely 
determined by the one-body correlation matrix Qjj — 
(c\cj) — Tr (gCjCj), where g is the density matrix. Now 
consider a bipartition of the full system into two subsys- 
tems A and B, and let i and j denote sites within A. 
Then G tj = (cJ Cj .) = (RGR T )ij where R u = 5 U is an 
oblong matrix of dimensions Na x (Na + Ng) with l's 
along the diagonal; Na(b) is the size of the A(B) sub- 
space. I.e. R spatially projects onto A. Thus G is a 
submatrix of G, and a key fact, due to PescheP, is that 
we may write Gj - = Tr (gc\cj), where g = exp(—K)/Z is 



the reduced density matrix (RDM) and K = c\ Cj is 
the dimensionless 'entanglement Hamiltonian' (both re- 
stricted to A). One then finds G = {exp(r T ) + 1} \ 
and the normalization Z follows from Tr(^) = 1. The 
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N A eigenvalues {j a } of r T = In (G _1 - l) are the en- 
tanglement 'quasienergies', and the eigenvalues of G are 
then Fermi functions of the quasienergies^, viz. 



fa 



exp(7 a ) + 1 



(4) 



For our system, the translation invariance along x means 
k x is a good quantum number, and for each k x , the sys- 
tem can be regarded as one-dimensional. Thus, Hu(k x ) 
and Qij{k x ) axe of dimension N y , and I and J label rings 
and not single sites. In our study, the Fermi energy Ef is 
always placed within some bulk gap. For periodic bound- 
ary conditions, Q is then a sum of projectors onto the oc- 
cupied bands. With open boundaries, there will be edge 
modes which cross the Fermi level. In either case, we take 
the A subsystem to be the bottom part of the cylinder, 
with y e [1, M]. Thus G is the upper left M x M block 

of g. 
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2.2.2. Rank of G and 1 - G 

The eigenvalues of G may contain exact zeros or ones. 
The number of zeros and ones, denoted as Do and D\ re- 
spectively, are by definition the dimensions of the kernels 
of G and 1 — G. If v bulk bands are occupied with pe- 
riodic boundary conditions, then the total rank of Q{k x ) 
is vNy/q, since each of the q bands contains an equal 
number of states. Thus if M > vN y /q, the rank of G 
will also be vN y /q. For M < vN y /q, the rank of G is M. 
Thus, 

rank(G) = min(M, vN y /q) , (5) 
D = M - rank(G) , (6) 

and similarly 

rank(l - G) = min^M,N y (l - v/q)\ , (7) 
D l = M - rank(l - G) , (8) 

where N y (l — vjq) is the rank of 1 — G, the projector onto 
unoccupied bands. 

As we shall see, with open boundary conditions, the 
rank of Q changes with k x whenever an edge state crosses 
the Fermi level. As a result, Do and/or D\ are discon- 
tinuous at such k x values for certain range of M. It is 
easy to verify that the condition for Do = Di = is 
M < mm(V/<7, 1 - vjq) X N y . 

2.2.3. Entanglement occupancy, quasienergy, and spectral 

flow 

Fig. [T] shows the energy spectrum and entanglement 
occupancy for p/q = 3/7, N y — 28, and M — 14, using 
open boundary conditions. The Fermi energy Ef lies 
within the third gap. The total Chern number of the 
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FIG. 3: (Color online) Hofstadter model with flux p/q — 3/7 
per plaquette on a cylinder of height N y — 70. Top: detail 
of energy spectrum showing lowest three energy bands. The 
Fermi level lies at Ef — —1.9 (grey horizontal line). The 
total Chern number of the occupied bands is C occ = 3. Red 
(blue) vertical lines indicate k x at which the lower (upper) 
edge modes cross the Fermi energy. Bottom: entanglement 
occupancies f a after tracing out the upper half of the cylinder. 



three occupied bands is C occ = X)j=i = 1; which is 
also the number of times the lower edge mode, localized 
at y = 1 (the red line in the plot), flows across the gap, 
i.e., the winding number^. The sign of C occ is reflected 
in the direction of the edge flow, e.g. C occ = — 1 for 4 
filled bands, and the lower edge flows downward. Several 
features are noteworthy: 

(a) Most levels are clustered near / = and / = 1. This 
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FIG. 4: (Color online) Same system as in Fig. [3] Top: entan- 
glement quasi-energies j a . 20 levels are below the large 'gap' 
- the same as the number of occupied levels in the full system. 
Colored levels are also plotted in Fig. [3] with the same color 
scheme. Bottom: quasienergy in 2D polar coordinates. The 
radius is the quasienergy and the polar angle is k x . The black 
circle corresponds to 7 = 0. There are three curves spiral- 
ing outward in the clockwise direction, corresponding to total 
Chern number C = 3. 



reflects the fact that Ef lies inside a bulk gap. Taking lin- 
ear combinations of the occupied states in the full system, 
one can create wavefunctions which are mostly confined 
to either (A or B) subsystem. The same consideration ap- 
plies to unoccupied states. In the thermodynamic limit, 
the fraction of occupied states in A and in B should be 
the same as that for the full system. This is confirmed 
by our numerical results. The entanglement eigenstates 



with / ~ i are localized along y « M. As is the case for 
the edge modes of the full Hamiltonian, the number of 
entanglement levels flowing between / ~ and / ~ 1 is 
the same as the total Chern number number C' occ of the 
filled bands. This is depicted in Fig. [T]for C occ = 1 and 
in Fig. [3(b)1for C occ = 3. 

(b) The occupancy f a (k x ) is discontinuous at k x = k x , 
where the lower open boundary edge modes of H cross 
the Fermi level; these are the red curves in Fig. 3(a) 



For example, in Figs. [T] and [2j where entanglement lev- 
els with different indices a are plotted in different col- 
ors, the tenth (light blue) and eleventh (magenta) level 
occupancies are each discontinuous at fc£, but satisfy 
/a=io(A£ + + ) = / Q =ii(fc^-0 + ). The number of distinct 
such k x values is the number of lower edge mode Fermi 
level crossings, which is |C occ |. Thus, as one increases 
k x through each such crossing, the number of levels with 
/ ~ 1, and hence the total occupancy J2 a fa of the A 
subsystem, drops discontinuously by unity due to the ex- 
clusion of the edge mode. Eventually the / ~ 1 levels 
are repopulated due to the aforementioned spectral flow. 
In order to conserve the rank of G{k x ) (in the case of 
M > rank(tj), cf. eqn.[6| upon increasing k x by 2-tt, then, 
there must be a discontinuous repopulation of the / ~ 
levels. This occurs when the upper open boundary edge 
modes cross E F ; these are the blue curves in Fig 



As these modes have a vanishingly small projection onto 
the A subsystem in the thermodynamic limit, they lead 
to no discontinuity in the total occupancy of A. A sim- 
ilar analysis of the occupancy discontinuity has recently 
been given by Alexandradinata et aly^* 

(c) From the occupancy spectrum, one can invert the 
Fermi distribution (eqn. [4]) to get the quasienergy spec- 
trum {7a}. A quasienergy plot more clearly reveals en- 
tanglement spectrum near / = and / = 1, where many 
levels are clustered. In Fig. |4(a)| (and equivalently the 
polar plot Fig. |4(b)[ ), some key features are apparent. 
First, a substantial number of levels are clustered at 
7 m 140 (/ ~ 10~ 60 ) and are separated from the re- 
maining levels by a pronounced gap. Actually this is a 
numerical artifact and these levels all lie at 7 = 00. Re- 
call the earlier result rank(G) = mm(M,i/N y /q) in eqn. 
[6] for a system with periodic boundary conditions. Here 
we have p/q = 3/7, N y = 70, M — 35, and v = 2 since 
Ef is placed in the gap between the second and third 
bulk bands. Thus we would expect rank(G) = 20, and 
since the row dimension of G is AI — 35, there should be 
15 levels with / = 0, corresponding to 7 = +00. Had one 
looked at a system with more than half filling, one would 
find entanglement quasienergies clustering at 7 = —00 
instead, where the entanglement occupancy is exactly 1. 
These would correspond to the kernel of 1 — G. 

(d) Since our system has open boundary conditions, 
there are edge states, and there is a discontinuity in the 
quasienergy spectrum at each k x and k x value. When 
both edge states lie below Ef, one counts 20 finite 
quasienergy levels. When k x lies between consecutive 
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FIG. 5: (Color online) Normalized eigenfunctions ijj a {y) of 
G for p/q = 3/7, N y = 112, M = 56, and v = 3 at k x = 
2ir/7 (open boundary conditions). The Chern number of the 
filled bands is C OC c = 1. Color corresponds to the sign of 
the wavefunction (red for positive, blue for negative), and 
intensity to amplitude (white for zero intensity). The black 
points are the entanglement energies y a . 



k x and k x values, one of the lower boundary edge states 
has crossed the Fermi level, and the rank of G decreases 
to 19. The spectral flow in the vicinity of 7 w is contin- 
uous, however. Discontinuities in the entanglement ener- 
gies occur for large values of I7I, where the occupancy is 
close to of 1. When an edge state passes from below Ef 
to above Ef, the rank of Q changes discontinuously by 
— 1. For increasing k x , this occurs at any of the three k x 
points in Fig. 



Such an edge state has almost per- 
fect projection onto the A subsystem, hence its depopu- 
lation leads to a sudden rearrangement of entanglement 
levels with large negative quasienergies j a (fa ~ 1) and 
a loss of one such level. For k x — where the change 
Arank(£/) = +1, the 'extra' level enters via a discon- 
tinuous rearrangement of the levels with large positive 
quasienergies (f a rj 0). (The situation is reversed if Ef 
lies within the first gap, in which case Arank(C?) = +1 
at each k x and Arank(Cf) = — 1 at each k x .) We see this 
clearly in Fig. 4(a) where the number of finite 7 levels 



changes from 20 to 19 when k x lies between consecutive 
fci and fc" values. 



2.2.4- Entanglement eigenfunctions 

In Fig. [5] we plot the eigenfunctions ijj a of G(k x — 
2tt/7) for a larger p/q = 3/7 system, with N y = 112, 
M = 56, and v — 3 at k x = 2t:/7. The Chern number 
of the filled bands is C occ = 1. States with finite j a 
are spatially resolved. The rank of G is vN y /q = 48, 



corresponding to states #9 through =#=56 in the plot. The 
dimension of the kernel of G is then dim(G) — rank(G) = 8. 
These states all have f a — 0, i.e. 7 a = +00, which is 
rendered as the flat ceiling of the black curve in the figure. 
They form the speckled region in the left of the figure. 
Note that states with 7 a ps (/ ~ \) are localized near 
the cut y = M, and that those with large 7 a are localized 
away from the cut. We shall return to this point later, 
toward the end of the paper, after we discuss Wannier 
center flows. 



2.2.5. Effect of changing M 

Since the magnetic unit cell is set along the y-direction, 
N y and M must be chosen as integer multiples of q if 
there are to be an integer number of unit cells in the full 
system and/or the lower (A) subsystem, respectively. As 
shown in Appendix [Aj changing N y to N y + m (m G Z) 
keeps the lower edge modes intact, but shifts the k x val- 
ues for the upper edge modes by —2irmp/q. For example^ 
the lower edge modes (red lines) are the same in Fig 
and [2j but the upper edge modes (blue lines) in Fig 
are shifted in k x by —6tt/7 relative to those in Fig.[l] 

It turns out that changing M affects the entanglement 
occupancy in the same way as changing N y would affect 
the edge modes. This is shown in the bottom rows of 
Figs, [l] and [2] We should mention that keeping M fixed 
while changing N y will not change the occupancy spec- 
trum in any appreciable way because that only shifts the 
upper edge modes. Changing N y will thus change the 
k x values, and consequently where the rearrangements 
of the / ss parts of the entanglement spectrum occur, 
but will not affect the spectral flow for 7 « 0. The reason 
will become more clear in the next section. 



3. ADIABATIC PUMPING OF BAND 
PROJECTORS 

The entanglement level occupancies f a are eigenval- 
ues of the restricted correlation matrix G = RQR T . In 
searching for an intuitive picture of the various features 
of the entanglement spectrum, it is then natural to ex- 
amine the unrestricted projector Q . We found that much 
information can be extracted from Q itself. 

In this section, we will use periodic boundary condition 
in the y-direction, i.e.. z — 1 and N y mod q = in cqn. [3] 
The Hamiltonian of eqn. [3] then satisfies 

H(k x +0) = TlH(k x )T y , (9) 

where <f> — 2irp/q as before and where T y is the transla- 
tion operator by one lattice spacing in the y-direction: 



T 



iV„-l 



(10) 



The unitarity of T y guarantees that the bulk bands re- 
peat themselves for q times over the interval k x £ [0, 27r]. 
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FIG. 6: (Color online) Full system band projectors for p/q = 3/7, N y 
k x — 2-Kn/q. Bj is the projector onto the j band, and Q„ = Bi + . . 
magnitude of the matrix elements are represented by intensity and their sign by color (red positive, blue negative, white zero). 
2/1 and y 2 are the row and column indices of the projectors. Blue rules mark boundary of the magnetic unit cells. Red rules 
mark the bipartite cut, so the top-left quadrant of Q corresponds to the restricted correlation matrix G. Only k = 1 and 3/2 
are shown here due to space restriction. Projectors at all other integer n can be obtained by shifting all matrix elements of 
those at k = 1 along the diagonal by t as the solution of eqn. [16] Similarly, projectors of all half odd-integer k are obtained by 
shifting those of k — 3/2. As k x is increased from to 2-7T (k from to q), the adiabatic pumping is evident in the diagonal 
motion of all matrix elements of the projectors. Note in particular that at integer k, each diagonal block of <J 3 consists of three 
sharply localized packets: the top-left one is contributed by B\, whereas the rest two result from the constructive addition 
of the diagonal blocks of B 2 and B3. Similarly, at half odd-integer k, Bi and B 2 add constructively, yielding the two sharp 
packets in each diagonal block of Q 2 . The Chern number corresponds to how many magnetic unit cell boundaries (blue and 
red rules) any diagonal matrix element has passed by in one pumping period. Equivalently, it is the sum of diagonal matrix 
elements which are transferred across any magnetic unit cell boundary in one pumping period. For the Q v type, it is intuitively 
how many packets are transferred. 



With each successive increase of k x by 2irp/q, the spec- 
trum repeats and the corresponding energy eigenstates 
are shifted by Ay = 1. Denoting Bj(k x ) as the projector 
onto the j th band, we have that Q(k x ) = Gv{k x ) is the 
projector onto the lowest v filled bands, 

V 

GM = J2 B oi k *) ■ (11) 

The covariance in k x and y is reflected as 

B 3 {k x + <f>; yi ,y 2 ) = Bj(k x ; y x + 1 , y 2 + 1) , (12) 

where y\ and 2/2 are row and column indices for Bj(k x ). 
Translational invariance on a scale of the magnetic unit 
cell corresponds to 

Bj(k x ; yi , y 2 ) = Bj(k x ; y t + q , y 2 + q) (13) 

The same relations hold for Q v . 



While these projectors are explicitly constructed us- 
ing the Bloch states, which are spatially extended, the 
fact that their eigenvalues are degenerate (either or 1) 
means one may construct localized eigenstates around 
the cylinder, for each k x , by recombining Bloch states of 
different k y with the same eigenvalue. In the continuum 
limit, where q — > 00 with p finite, these correspond to 
the familiar Landau strip basis. In fact, the projectors 
themselves are localized: An illustration is provided in 
Fig. [6j which shows several Bj and Q v for p/q = 3/7 at 
k x = 2nK/q for n = 1 and k = 3/2, both of which are 
local extrema of the energy bands. The size of the mag- 
netic unit cell naturally divides the projectors into blocks 
of size q x q. It is not surprising that the off-diagonal 
blocks drop exponentially, a consequence of the analytic- 
ity of Bj in complex kjp**. What is perhaps unexpected 
is that at band troughs (integer k for odd j and half-odd- 
integer n for even j in Fig. [61), the projectors Bj, and es- 
pecially Q v , are quite well- localized even within the diag- 



onal blocks. The diagonal matrix elements of the projec- 
tors correspond to electron density at the corresponding 
y coordinate. For single bands (Bj), any q consecutive 
diagonal elements sum to 1, thus one may think of them 
as constituting a wavepacket, and the projector Bj as 
consisting of N y /q such wavepackets (one per magnetic 
unit cell). For each wavepacket, the weight is dominated 
by one or two elements, as one can see in Fig. [6j The 
localization of the projector sums Q v is even more promi- 
nent: when the gap between two neighboring bands is at 
a minimum, their projectors add constructively, resulting 
in two sharply localized dots on the diagonal line of Q v . 
For example, in Fig. |6j B 2 and B 3 add constructively at 
k — 1, yielding the lower two dots in each diagonal block 
of C/3, and similarly, B\ and B 2 add up to Q 2 for k = 3/2. 
Now, q consecutive diagonal elements in Q u must sum to 
v, thus each of the v dots can be intuitively understood 
as one localized wavepacket. The constructive superposi- 
tion of neighboring bands then indicates the correspond- 
ing single-band wavepackets have opposite parity so that 
the off-diagonal elements cancel each other. To relate to 
the aforementioned "strip" states, we note that any col- 
umn of a projector is an eigenstate of the same projector, 
with eigenvalue 1^. The diagonal nature of these band 
projectors thus ensures the existence of such strip states. 

We now establish a connection between the band pro- 
jectors and the seminal work of Thouless et al\^ on the 
Chern numbers for the Hofstadter bands. Consider first 
the individual band projectors Bj. We write 



y = q£ + m 



2ttk 

q 



(14) 



where £ and m are integers. Thus £ is the magnetic unit 
cell coordinate, and m the coordinate within each such 
cell. For a single band, denote the position of any of its 
wavepackets as m(«;), then (k x ,y) covariance implies 



m(/i + tp) = m(/t) — t 



t G Z 



(15) 



Of course both m and n are only defined modulo q. The 
relevant quantity in the k x pumping is the 'velocity' of 
the packet (with k x as 'time'), i.e., the number of sites it 
traverses when n is effectively increased by 1: 



tp = sq + 1 , s € Z , |t| < q 



(16) 



A graphical construction is shown in Fig. [7j Clearly, t 
will be the number of packets transported through any 
given boundary during the cycle k x — » k x + 2tt (number 
of blue flow lines in the figure). It is also equal to the 
number of magnetic cells traversed by a single packet. 
There is however a mod q ambiguity associated with the 
sign indeterminancy of t, e.g. for p/q = 3/7, one has 
that (t, s) = (5, 2) and (- 



16 In 



-2, —1) both satisfy eqn. 
general, without looking at intermediate k snapshots, one 
cannot tell if the packet had advanced by t or retreated 
by q — t. We have examined different p/q ratios on both 
square and triangular lattices, and we find that for the 
lowest band on a square lattice, the ambiguity can always 
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FIG. 7: (Color online) (m, k) construction for p/q — 3/7, 
assuming m = when k = 0. Red line: flow in the order 
of T y translation whereby k — s> k + p and y — > y — 1. Blue 
lines: flow in the order of k x pumping. The number of time 
a packet is transferred across any boundary line, e.g. m = 
line, is the same as the number of blue flow lines crossing the 
boundary line. There is a mod q ambiguity as can be seen 
from the validity of both the solid and dashed blue flows. 



be resolved, without needing to inspect intermediate k, 
by picking the value of t which has the smaller magnitude 
\t\, e.g. t = —2 instead of 5 for B\ in Fig. [6] Intuitively, 
this means the packet moves toward the nearest possible 
position allowed by eqn. [15] 

Eqn. fl6] is recognized as the Diophantine equation of 
TKNI 1 ^ with r = 1, according to which t is simply Ci, 
the Chern number of the lowest band. The heuristic of 
taking the smaller \C\ \ in resolving the mod q ambiguity 
agrees with ref . [T] Since this picture does not distinguish 
between different bands, the Chern numbers of all bands 
are equivalent mod q. 

The correlation matrix Q v for v occupied bands has v 
packets per diagonal block, with each moving according 
to eqn. 16 as required by the (k x ,y) covariance. How- 
ever, their collective motion depends on their relative 
spacing. Consider for example the diagonal blocks of 
G3 for k — 1, as shown in top right panel of Fig. [6] In 
each diagonal block, the diagonal elements with domi- 
nant weight (packets) are at y = (2,4,7). Then (k x ,y) 
covariance requires that at k = 2, they are moved to 
y = (2 - t, 4 - t, 7 - t) = (4, 6, 9) where t = -2 is the 
solution to eqn. [16] as discussed before. Since y = 9 is 
simply the y = 2 element of the next diagonal block, 
thus in each diagonal block, the packets of k = 2 are 
at y — (4,6,9 mod 7 = 2). Comparing this with those 
of K = 1, one can see that, effectively, only one packet 
moved from y = 7 to y = 6. We found that for any 
v 7^ q, this observation holds true (for v = q, Q v is iden- 
tity). That is to say, as n —> k + 1, the net effect is for 
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only one of the v packets to change position. We will not 
attempt to explain this observation, but rather take it as 
a starting point, and explore its implications. Mathemat- 
ically, this observation - the reduction of the motion of 
multiple wavepackets to that of a single mobile packet at 
a time - means that the positions of dominant diagonal 
elements (the packets) can be labeled in such a way that 



m,-(K + 1) = m i+ i(/c) 



i=l,2, 



,v-l 



(17) 



with each r/ij still satisfying eqn. 15 
mj(/t) — t. 



i.e. 



mi{n + tp) = 

Note that this does not mean rrii can be iden- 
tified with the wavepacket of a single band: they result 
from constructive superposition of single-band projec- 
tors, as discussed earlier. To illustrate eqn.[T7| take again 
G3 as an example: at K = 1, (mi, TO2, 1713) = (7,2,4), 
while at k — 2, (mi, ma, 7713) = (9 mod 7 = 2,4,6). 
Then as k increases by 1, the effective change is of one 
packet (the mobile one) moving from mi(«) to m„(ft + 1) 
with stride t v , 

m„(/t + 1) = mi(/t) - t v , \t v \<q. (18) 

The RHS is therefore mi (ft + t„p), while the LHS is 

m^K + 1) = m u -i(K + 2) = • • • = mi(K + 1/) (19) 

from eqn. |17| Thus t v is determined by 



hence a symmetric flow from to 1 with its wavefunction 
localized at the opposite end. (Entanglement occupancy 
with periodic y boundary is shown in Fig. [9] in the next 
section). Furthermore, If we change the position of the 
entanglement cut M (not necessarily along a magnetic 
cell boundary, for example), this will simply change the 
value of k x when a packet hits the cut, whence the k x 
translation shown in Figs. [T] and [2] 

While the entanglement spectrum only reveals the to- 
tal Chern number, the correlation matrix retains some 
information about the individual Chern numbers of con- 
stituent bands, manifested as the separation between its 
wavepackets. Note that eqn. 18 can be taken as a defi- 
nition of t v with arbitrary v < q, without interpreting v 
as the number of filled bands. After all, the total Chern 
number of the lowest v bands is the same whether or not 
they are filled. We explicitly replace v with n < q below 
to avoid any such connotation. From eqn. |17| and |18| we 
have 



l n+l 



(k) - m n (n) = t, 



t, 



(21) 



t v p = sq + u 



s G Z , \t v \<q 



(20) 



Again, there is a mod q ambiguity because of the sign 
indeterminancy of t v . For square lattice, the heuristic of 
using the smaller \t v \ still seems to hold, e.g. while both 
(*3,s) = (1,0) and (*3,s) = (—6,-3) satisfy eqn. 20 the 
actual system picks t^ = 1. Ref. |2T)1 mentioned that s and 
t v cannot simultaneously be odd for cither the hexagonal 
or triangular lattices. Incidentally, for v = p, t p = 1 is 
always a solution with the corresponding s — 0, i.e. the 
total C her n number of the lowest p bands is always 1. 

Eqn. 20 is the TKNN Diophantine equatiorfQ for r = v. 
There, t v is the total Hall conductivity (the sum of the 
Chern numbers) of the v occupied bands. It is also the 
winding number of the energy edge states in the u th garM 
We now have a third interpretation: it is the number of 
sites traversed by the mobile packet during each k incre- 
ment. Equivalently, it is the number of mobile packets 
transported across any magnetic unit cell boundary dur- 
ing the cycle k x — > k x + 2n. 

The entanglement spectrum can now be understood 
intuitively. Whenever the mobile packet leaves the lower 
half-cylinder through the cut between M and M + 1 (in 
Fig. [6j proceeding from top-left quadrant through the red 
line into the lower-right quadrant), there is an occupancy 
flow from / = 1 to 0. The number of flow lines is then 
equal to the number of packets which move through the 
cut, which is the total Chern number. In the periodic 
y boundary case, where the cylinder is compactified into 
a torus, the flow across M = N y /2, is always concomi- 



thus the two packets at m n +i and m n are separated by 
a distance of C n . Quantities such as the four-point cor- 
relation IF(A) = (c y c y c y+A c y+A ) thus have peaks at 
A = Cj apart from A = q, 2q, etc. 

On a square lattice, the Hofstadter model exhibits a 
particle-hole symmetry. This implies that Cj = C 9 +i_j. 
For eve n q, t he bulk spectrum is known to have no cen- 
tral gap 1 21 * I, therefore the Chern numbers of the two 
central bands are not individually well defined, and one 
can speak only of a Chern number for the pair. It is in- 
teresting to notice its implication on the distribution of 
the wavepackets within each unit cell: if on the contrary 
there is a central gap, then ti — C\ + C2 + ■ ■ ■ + Ca — 
C q + Cq-i + - ■ ■ C|+i' Since the total Chern number of all 

bands must be zero, we must have ti = 0. Now accord- 

2 

l+i( K ) 



ing to eqn. fl7l and 
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mi(re) — t± = 7?t,i(k), 
so the (| + l) th packet and the first one are forced onto 
the same site. Thus, the fact that there is no central gap 
in this case guarantees that there will be no packet 'colli- 
sions'. The natural question to ask next is how the wave 
packets are arranged when v > q/2. To illustrate this, 
consider a specific case with p/q = 5/8 and v = 7 filled 
bands. The Chern numbers of the lowest three bands are 
—3, 5, and —3, respectively, and particle-hole symmetry 
guarantees that these values repeat for the upper three 
bands. The central two bands therefore have a combined 
Chern number of C4 5 = 2. Since Chern numbers rep- 
resent the separation between wave packets, we can fill 
in the first four packets with no difficulty (the position 
of the first packet being arbitrary) . The location of the 
fifth packet cannot be determined because C4 is not well- 
defined, but the location of the sixth packet is found by 
shifting the fourth one by C4 5 = 2. The rest of the pack- 
ets can be filled in a similar fashion. Thus the vanishing 
of the central gap implies an indeterminacy of the po- 
sition of the (§ + 1) wavepacket. To resolve this, one 
tant with another packet moving from y = N y to y = 1, can add in an infinitesimal second-neighbor hopping that 
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breaks the particle-hole symmetry and results in a small 
central gap. For example, one can introduce a second- 
neighbor hopping t' along one of the two diagonals in 
each unit cell (say in the direction x — y). This construc- 
tion interpolates between the square lattice when if = 0, 
and the triangular lattice when t' = 1 (see Appendix |b|). 
For t' <C 1, we find that the Chern numbers of the central 
two bands are resolved as C4 = 5 and C5 = —3. C4 can 
now be used to determine the position of the 5 th packet. 



WANNIER CENTER FLOW 



The Diophantine equation (20) describes a mod-g 



property of the Hofstadter problem, which is a result 
of the (k x ,y) covariance of eqn. 12 No knowledge of 
intermediate k x € [k</>, (k + l)<p) is required in obtain- 
ing eqn. |20| This comes at a price of the ambiguity in 
t (mod q) and s (mod p), which intuitively contain the 
information of the direction in which any given packet 
is moving. To settle that, we need localized eigenstates 
of the projectors over the full range of k x , namely the 
Wannier functions. 



4.1. Wannier functions in ID 

The application of Wannier functions to the analysis of 
topological band structures has recently been developed 
in refs. Q31 [THl an d HS1 Consider a ID periodic system 
with N unit cells and q bands. The energy eigenstates 
are 



where A k — {ipk\ id k \ipk)i an d where A is quantized ac- 
cording to 



2tt 



1 



Aj = -- ldkA k +I , JeZ . (26) 
2ir 



Thus for a single band, the state \&\ r ) is localized at unit 

cell I with an offset 7/271", where 7 = J dkA k . 



When the internal space of \ip n . k ) coincides with real 
space {e.g. the lattice site m within the magnetic unit 
cell in Hofstadter problem), one may also introduce the 
cell functions \u n k ), 



K fc ) = v k , v k = J2 M z~ lhmlq H ( 27 ) 



such that 



(x;m\*n,k) =e ik < x+m M (m\u n , k ) 



(28) 



In this case, an alternative definition of the position op- 
erator is 



R x = X © — , rh 

q 



V] m\m)(m\ , (29) 



where rh measures the position within the unit cell. This 
amounts to replacing \ip n ,k) with \u n ^) in the Wannier 
center calculation, 



I* 



n.k) 



\k) ® |i/Vfc) 



(22) 



where (x\k) = e X /\N is the Bloch phase of unit cell 
x, n is the band index, and \ip n ,k) is the n th eigenstate 
of the q x q Fourier-transformed Hamiltonian H(k). The 
Wannier functions can be defined as the eigenstates of 
the projected position operator^ PRP, where P is the 
projector onto one or several bands, and R = X (g) 1 is 
the position operator. The eigenstates are linear combi- 
nations 



|4> 



(23) 



satisfying PRP\<S> X ) = A |$ A ), with /fe +27r = f k . In the 
bulk, one may use the substitution X — > — idk on the 
{|fc)} states. It is not valid near the periodic boundary 
x = 1 and N . Upon an 'integration by parts', one obtains 



(A - A k ) f k = i 



dfk 
dk 



(24) 



The solution is 



7n -> In = / dk (u n , k \id k \u n> k) 



(30) 



= Jn + q 1 / dk (i/j n ,k\rh\ip n , k ) . (31) 


The extra term is the average position of the band within 
the unit cell. 

Incidentally, by dropping the Berry connection in 
eqn. |25[ one obtains the conventional Wannier function, 



v> = ^-'"l*. 



Using X — i dk, it is easy to check that 
( aj |A-®l| 0j ) = ~+/ 

Z7T 



ftl 

( ai \x®-\ ai ) 



2tt 



(32) 

(33) 
(34) 



fk = fo exp 



dk\A k , - A) 



(25) 



i.e., it has the same spatial center as |A/) (but they differ 
in their spatial spread). This is in agreement with the 
result of ref. [251 
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For multiple bands, the eigenvalue problem becomes 
|*a}= ' PRP\^x) = X\^x) (35) 

n,k 

where n runs over v involved bands, again subject to 



periodicity j n ,k+2n = fn- Similar to eqns. 24 and 25 the 
solution ia^S 



fn,k — 6 %Xk Wnn' fn',k=0 , 



(36) 



where 



W = P exp 



(37) 



i J dk'A(k') 
. o 

is the Wilson loop, P denotes path ordering, and where 
A mn (k) = (tp m ,k\idk \ipn,k) ■ (38) 

is the Berry connection (geometric vector potential). 
Since A is Hermitian, W is unitary and its eigenvalues 
may be written as e i9 ™ . The eigenvalues of PRP are 



A,, 



2tt 



I 



(39) 



For two-dimensional lattices, Wannier functions can be 
defined at each k x . For a single band n, its y center is 
the Berry phase j n (k x )/2ir. The band Chern number is 
the negative of the winding number of "f n (k x ) over the 
interval k x £ [0, 2ir] : 



C?7 



1 

2~7ri 



2m 



d\V k x (ip„(k)\ V fe |V>n(fc)) 



dk v (ip n (k)\d ky \ipn(k)} 



k x =0 



7n (0) - 7n (27r) 
2ir 



(40) 



Thus the Wannier center shifts by — C n (magnetic) unit 
cells over k x — > k x + 2n, as found by Qi in ref. H3J We 
have seen in Sj3jthat the packet associated with band n is 
translated by — C„ lattice sites during k x — > k x + 2ir/q, so 
over k x — > k x + 2ir, it will be translated by —C n magnetic 
unit cells, in agreement with the Wannier picture. 
For multiple bands, we have 



hence 



2tt 

del U = exp {i j dk I r . i,. 



°™ = ^ 

kj— 1 n— 1 



(41) 



(42) 



Then similar to the single band case, we conclude that in 
2D, the total shift of all (inequivalent) Wannier centers 
is given by the sum of the Chern numbers. This is re- 
flected in f|3]as the number of mobile packets transported 
through any given magnetic unit cell boundary. 

When there is no level crossing among {\ w .i} over the 
period of k x , one can combine the w and / indices. De- 
fine a composite index (i(w, I) = w+vl, with X w j — > A M . 
Sending k x from to 2tt amounts to an index shift, which 
is universal for all 9^ since there is no level crossing. Then 
the eigenfunctions 9 Al (fc a; ) of W (here taken to be peri- 
odic in the index modulo q) satisfy 



m (2tt) = e M+(T (o) 



cr e z 



(43) 



which is just a cyclic permutation in the w index with an 
offset a. One may think of the set of {e iS ™ } as v points 
on the unit circle where different indices / are equivalent. 
Then during the a cyclic permutation, the perimeter of 
the circle is covered by these v points for a times, i.e., 



± £ [<V(2tt) - e wJ (0)] = a 



(44) 



w— 1 



But according to eqn. 42 the LHS is simply the total 
Chern number of constituent bands: 



(45) 



4.2. Wannier center flow in Hofstadter problem 
and a general relation with the entanglement 
spectrum 

A slightly different definition of Wannier functions is 
used in our study of the Hofstadter problem to accommo- 
date the use of finite N y and periodic boundary condition 
in y. We follow here the recent work of Yu et al in ref. 
IT51 Instead of PRP, we numerically diagonalize 



(46) 



and use its eigenstates with non-zero eigenvalues as Wan- 
nier functions. States with zero eigenvalues are those 
projected out by P. For a single band, P = B n , while 
P = g v for v filled bands. e i2Tr */^ 
periodic-boundary version of 1 — i 



is essentially the 
dk and plays the role 
of — i dk in §4.1| An eigenvalue problem can be similarly 
solved, whose N — > oo limit is the same as those ob- 
tained before. Notice also that y corresponds to X © y 
as in e qn. [29] Assuming the eigenstates and eigenvalues 
of eqn. |46| are {\9)} and {e lS }, we compute three slightly 
different Wannier centers, 



Ve 



-N„ 



2ir 

(y) e = (9\ij\9) 
(Y) B = {6\Y®1\ 



(47) 
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k = k x q/2ir 



(b) <y> 

FIG. 8: (Color online) Wannier levels of the lowest band for 
p/q — 3/7 with 4 magnetic unit cells, (a): Using lattice co- 
ordinate y as position operator. Colored dots: yg. Notice 
that points of the same type/color are related via {k x ,y) co- 
variance, in this case advancing byp = 3inKasi/->i/ + 6 
(equivalent to y — 1). Colored lines: {y)e- Black-white back- 
ground: diagonal elements of Bi, black = 1, white = 0. See 
also Fig. [6] (b) : (Y) g . Corresponding levels have the same 
color as in (a). 
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(a) 2 filled bands 
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FIG. 9: (Color online) Wannier centers vs. entanglement oc- 
cupancy for p/q = 3/7 with different filling fraction and 4 
magnetic unit cells in the full system and 2 in the half sys- 
tem. Top: Wannier centers using lattice y coordinate, sim- 
ilar to Fig. 
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Each type/color of point corresponds to a 
mobile packet (because they are not related by the (k x , y) co- 
variance). In (a), advancing k by p = 3 does not change y by 
— 1 or 6, and in (b), advancing n by 1 does not change y by 2 
or —5. It is also clear that within each magnetic cell, only one 
mobile packet moves as/t-> k+1. with k = half-odd-integer 
in (a) and integer in (b). Center: W annier centers using mag- 
netic cell coordinate, similar to Fig. 8(b) Bottom: entangle- 



ment occupancy. Color indicates level index and is not related 
to those used for the Wannier centers. Notice the similarity 
between the Wannier center flow and the entanglement flow: 
a mobile packet crossing the magnetic cell within the bulk has 
almost the same shape as the entanglement downflow; while 
crossing from y = N y to y — 1 has the almost the same shape 
as the entanglement upflow. The plateau-like feature in the 
entanglement flow lines can be traced back to the y plot in the 
top row as the motion of a mobile packet within one magnetic 
cell. 



where 9 is wrapped in such a way that y s is restricted 
0.5]. 



to \0.5,N V 



Fig. 



shows the Wannier centers 
defined above for the lowest band of p/q — 3/7 with 
N y = 28. y g (colored dots in (a)) has the proper trans- 
lational property: the dots are generated by shifting any 
single flow by successive multiples of q. In the vicinity of 



half-odd-integer k, a Wannier center migrates from one 
site to the next one which is —C\ = 2 sites ahead. In 
terms of the corresponding wavefunction (not plotted), 
what happens is that around integer k, it has a single 
peak at the site given by its (rounded) eigenvalue ye- As 
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FIG. 10: (Color online) Normalized eigenfunctions ip a (y) of 
QRQ for p/q = 3/7, N y = 112, M = 56, and v = 3 at 
fca; = 27r/7 (open boundary conditions). The Chern number 
of the filled bands is C OC c = 1. Color corresponds to the sign 
of the wavefunction (red for positive, blue for negative), and 
intensity to amplitude (white for zero intensity). 



k slowly moves toward the next half-odd-integer value, 
some weight is transferred to the next site, causing the 
eigenvalue yg to interpolate between the two sites. As k x 
is increased by 2n, q such migrations are made, i.e., each 
Wannier center is shifted backwards by C\ magnetic unit 
cells. In the bulk, (y) g (colored lines in (a)) and y g over- 
lap. This corresponds to the regime where X — > — i dk is 



valid in deriving eqn. 24 Near y — N y , part of the weight 
is pushed over the boundary to the y — 1 end, thus {y) g 

starts to deviate from y g and drops, until all weight is 
transferred to the other boundary. If instead of eqn. |46j 
the eigenstates of PyP are used, then the behavior in 
the bulk is still the same, but near the edge there will 
be avoided crossings in the Wannier center flow of (y) g . 
The flow of the magnetic unit cell coordinate (Y) g (col- 
ored lines in (b)) is similar to that of (y) g , but with an 

emphasis on the occasions when (y) g crosses a magnetic 
cell boundary. 

The Wannier centers of two and three filled bands are 
shown in Fig. [9] As k increases by 1, only one Wannier 
center per magnetic unit cell flows by a = — Yli=i Q 
sites (Eq. 45]), corresponding to the motion of a mobile 
packet in ^3p<7 = —3 for two filled bands and 1 for three 
filled bands. There is no level crossing in y g , hence eqn. |45| 
holds true: following any flow line from k x = to 2-7T 
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FIG. 11: (Color online) Normalized eigenfunctions ip a {y) of 
RQR for the parameters given in the caption to Fig. |10| 
The entanglement energies are plotted in the top half. Grey 
points correspond to j a — 00 (/„ =0). As expected, states 
with entanglement energy 7 a « are localized in the vicinity 
of the cut. 



leads one to —a levels beneath the starting point. Notice 
the same behavior in the entanglement quasienergy of 

Fig.gr 



There is a striking similarity between the magnetic 
cell coordinate flow (Y) and the entanglement occupancy 
spectrum, shown in the center and bottom rows of Fig. [9] 
the upward flow of (Y) looks exactly like the upward 
flow of /, while the downward flow of (Y) within each 

Y —> Y — 1 sector looks like the downward flow of /. This 
can be understood in the following way: the spectrum of 

Y can be thought of as a coarse-grained version of y, and 
should look just like (Y) and (y) shown in Fig. [9j so the 
effect of coarse-graining is to suppress the flow within a 
coarse-grained cell, and enhance the flow migrating be- 
tween different cells. One may push the coarse graining 
to the extreme where all sites with y < M count as y — 0, 
and all above M as y — 1. Then the only significant flow 
is from y = to 1, corresponding to Y = 2 to Y = 3 in 
Fig. [9j and from y = 1 to y = 0, corresponding to Y = 4 
through the periodic boundary to Y = 1. In Fig. [9] such 
a coarse graining would keep the upward Y flow intact, 
while push Y = 3 and Y = 2 lines to top and bottom 
respectively for the downward flow, making it look just 
like the entanglement occupancy spectrum. In fact, one 
can prove that the coarse-grained Wannier spectrum is 
identical to the entanglement spectrum: Consider two 
arbitrary projectors P and R. One can think of them 
as two matrices of the same dimension (zero-padded, if 
necessary, to fill out the dimensions). We claim that 
PRP and RPR have identical eigenspectra. To see this, 
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assume \%jj a ) is an eigenstate of PRP with non-zero eigen- 
value A„. Then 



\1> a ) = ^-PRP\i> a ) 



(48) 



and therefore P\tp a ) = \ipa)- Thus, \ip a ) is an eigenstate 
of P with eigenvalue 1. It then follows that 



RPR |Va) = KR IV'a) 



(49) 



from which it follows (using R 2 = R) that \ip a ) = R\if> a ) 
is an eigenstate of RPR with the same eigenvalue A a . 
Thus the non-zero spectrum of PRP belongs in that of 
RPR, and vice versa, so they are identical. Since the 
coarse-grained position operator R is (the complement 
of) the projector used in constructing the restricted cor- 
relation matrix G, the entanglement spectrum of G = 
RQR is identical to the coarse-grained Wannier centers 
QRQ, and the entanglement eigenstates are obtained by 
projecting the coarse-grained Wannier states onto the rel- 
evant half space. 

In Fig. 10 we plot the normalized cigenfunctions ^a{y) 
of QRQ for the case p/q = 3/7, N y = 112, M = 56, and 
v = 3 at k x — 27r/7. Note how the behavior of the 
eigenfunctions mimics that of the Wannier states, with 
ip a (y) localized at a position which moves across the en- 
tire cylinder as the label a advances through a range cor- 
responding to the rank of QRQ (f x 112 = 48 in this case; 
see £ 2.2.2). States #1 through #64 belong to the kernel 



of QRQ and are all degenerate. The insertion of the real 
space projector R thus fails to resolve these wavefunc- 
tions in real space, which explains the speckled pattern 
on the left side of the figure. The Wannier states are 
better localized however, since R may be considered a 
coarse-gr aine d approximation to y. 



In Fig. 
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we plot the normalized eigenfunctions ip a (y) 
of RQR for the same parameters, along with the corre- 
sponding entanglement energies j a - (If we remove the 
upper half of the cylinder, where the wavefunctions van- 
ish, this is a repeat of Fig. [5j) According to our defini- 
tions, 



\i> a ) = R\i> a ) I , v / WlW 



(50) 



Note how states of large positive j a (f a « 0) as well as 
states of large negative 7 Q (f a ss 1) are localized far from 
the A/B boundary. 



5. SUMMARY 

We have studied the entanglement spectrum and Wan- 
nier center flows of the Hofstadter problem. Most of the 
data presented in this paper was for the square lattice 
with p/q = 3/7 flux quanta per unit cell, but most of 
our observations are robust with respect to changing lat- 
tices, fluxes, and fillings. The entanglement spectrum of 
a subsystem exhibits spectral flow similar to that of the 



full system's energy edge modes: the total Chern number 
controls the number of flow lines, and its sign tells the 
direction of the flow. When open boundary condition is 
used in the full system, the entanglement spectrum ex- 
hibits level index discontinuity on the flow line. This is 
a manifestation of the crossing of the Fermi energy with 
the full system edge modes, which results in a total oc- 
cupancy discontinuity. Changing the location of the en- 
tanglement cut shifts the entanglement spectrum. This 
reflects the k x y covariance of the Hamiltonian: changing 
k x to k x + 2np/q is equivalent to shifting the system in 
y by Ay = 1. 

The behavior of the entanglement spectrum can be un- 
derstood by looking at the full system band projectors. 
These projectors are well localized and thus represented 
by wavepackets on their diagonals. For single bands, the 
packets flow under k x pumping. The (k x ,y) covariance 
then imposes restrictions on possible flow rate, described 
by a Diophantine equation first derived by TKNInP. Since 
the Chern numbers are also given by the same equation, 
the topological property of the system can be described 
equivalently in terms of the motion of these packets: the 
number of magnetic unit cells traversed by each packet 
during one period of k x is given by the Chern number, 
with the direction given by its sign. For multiple bands, 
the flow is that of the mobile packets moving under k x 
pumping, and the number of mobile packets crossing a 
given boundary gives the total Chern number of filled 
bands. The entanglement spectrum can then be under- 
stood as a measure of detecting when these packets cross 
a particular boundary, namely the entanglement cut. 

Using the (k x , y) covariance alone (and hence the Dio- 
phantine equations) only fixes the flow and the Chern 
numbers up to mod q because it only relates different 
k x points of fixed separation of 2ir/q. The localization 
of the projectors suggests the use of Wannier functions 
for smooth interpolation between these k x points. For 
single bands, the Wannier center at each k x is given by 
the corresponding Berry phase and is represented by one 
packet in the projector diagonal. The flow of the Wan- 
nier center is then described by the winding number of 
this Berry phase, which is the band Chern number. For 
multiple bands, the Berry phase is replaced by a set of 
eigenvalues of the Wilson loop operator. If there is no 
level crossing over the full range of k x , then all levels ex- 
perience a universal index bump of a as k x — > k x + 2tt, 
and a is given by the sum of Chern numbers. In com- 
puting the Wannier center, the position operator can be 
either the (magnetic) unit cell coordinate alone, or one 
that also includes the internal coordinates (lattice cell 
within each magnetic cell). The spectrum of the former 
is a coarse-grained version of the latter. One can take the 
coarse graining to the extreme of a bipartition, at which 
point the position operator becomes a real space projec- 
tor, and the coarse-grained Wannier spectrum becomes 
identical to the entanglement occupancy spectrum. 
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Appendix A: Incommensurate edge spectrum 

Here we first briefly review the edge spectrum with 
commensurate N y and q as studied in ref. |31 and then 
extend its argument to the incommensurate case. 



1. Review of commensurate edge spectrum 

The Schrodinger equation corresponding to the matrix 
equation H(k x ) yy >ip y i — eip y is 

-1py-l - lpy+1 - 2 COs(fca; + y4>) Ipy = ^Py (Al) 

cast into transfer matrix form, we have 



ipy 



y+1 = M„ 



tpy J \^y-l 

e — 2 cos(k x + y<fi) — 1 
1 



(A2) 
(A3) 



Notice that M y depends on e. The following boundary 
condition is required for eqn. |A2| to also cover the cases 
of y = 1 and N y , 



% = Vx +i = o 



(A4) 



It is then easy to verify that 

^g+i = [Gnf V>1, ^,=0 (A10) 

where I — 1, 2, . . . , L, hence the solution is an edge state 
exponentially localized at y = N y if \Qn\ > 1, and at 

y = 1 if I Qui < 1. 

The edge spectrum {e} satisfying the condition 
Q2i(s) — is the same as the full spectrum of a 
(q — 1) X (q — 1) system, so numerically the edge spec- 
trum of H(k x ,N y = Lq — l,z = 0) can be solved by 
diagonalizing its upper-left (q — 1) x (q — 1) submatrix. 

Note that eqn. |A10| implies the edge states, with ipo 
included, has a direct product form 

fQ°u\ 



Gii 
VGiTV 



ih 

\lpo-J 



(All) 



Then 



i.e., iptq+m = S'l^m with I = 0, 1, • • • L — 1 and m — 
0, 1, • • • q — 1. The iV-component magnetic cell part dic- 
tates the real-space behavior. In this case it is exponen- 
tially localized at either end. The g-component internal 
part is obtained by prepending ifjQ — to the solutions of 
the (q — 1) x (q — 1) upper-left block of H . This is by no 
means a general form of edge states, but we do also no- 
tice a similar decomposition in the zigzag edge modes of 
the Haldane modePSl. Note also that all Bloch states 
have such a decomposition, \^!{k,n)) = \k) <g> \tp n (k)} 
where (y\k) = e lky /vN is the Bloch phase and \ipn(k)) is 
the n th band eigenstate of the Fourier transformed q x q 
Hamiltonian. One may then say that — ilog(Qn) is the 
imaginary Bloch vector, and which of the UHP or LHP 
it resides in tells the localization of the edge states. 



2. Incommensurate edge spectrum 



<A5) 

My = M y M y _x ■ ■ ■ Mi (A6) 

and eqn. A~4 implies Afjv H +i is a triangular matrix, 

[Mn v+ i] 21 =0 (A7) 

The spectrum {e} consists of all energies satisfying 
eqn. \K7\ 

Notice that M y+q = M y , so when N y + 1 = qL with 
integer L ("commensurate"), 



M Ny+1 = Q L , Q = M q 



(A8) 



Now, products of up-triangular matrices are still up- 
triangular, so eqn. |A7| is satisfied if Q is up-triangular, 



Q 2 



(A9) 



In the thermodynamic limit where N y — > oo, one can 
extend the commensurate argument to incommensurate 
cases, Ny + 1 = Lq + m, with m = 0, 1, 2, . . . , q — 1. 

First, we note two properties of the transfer matrix, 

det(M v ) = 1 (A12) 
M y+m (k x ,e) = M y (k x +mcj),e) (A13) 

both are straightforward from definition. Eqn. |A13| 
expresses the same k x y covariance as eqn. 15 The 
(Lq + m)-step transfer matrix can then be divided in 
two ways, 

M Lq-\-m 

(k x )=M m (k x )Q L (k x ) (A14) 

= Q L (k x +m4>)M m (k x ) (A15) 

If Q(k x ) satisfies the commensurate edge condition 
eqn. |A9[ then 

Q L (kx) = ([Qn\ L ^^ L ) f (A16) 



[Q 2 
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where x is some number of no interest. We then have 



M Lq+m{k x ) 



An A12 

A 2 l A 2 2 



[Qn\ 




X 

[Q22] 



hence 



[M Lq + m } 2 i (k x ) = A 21 [Qnf (k x ) 



(A17) 



(A18) 



= A 2 i [Q 22 ]- L (k x ) 



where the second equality follows from det Q — 1, a con- 
sequence of eqn. |A12| 

Similarly, if Q{k x + mifi) satisfies the commensurate 
edge condition eqn. |A9[ we have instead 



M L q+m(k x ) 



Qn] 




[fia 



k x +m<f> 



An A i2 
A 2 i A 22 



(A19) 



thus 



[M Lq+m ] 2 i (M = [Q22] {kx + m<t>)A 21 (A20) 
= [Qn]~ L {k x + m4>)A 2 i 

We can then conclude that in the limit L — > 00, 

1. If {k x ,e) is a solution of the commensurate case 
N y + 1 = Lq at the lower edge, [Qn] {k x ) — > 0, 
then by eqn. A18[ it is also a solution of incom- 
mensurate N y + 1 = Lq + m. For y £ [l,Lq — 1], 
the wavefunction ip y coincides with that of the 
commensurate case, and in the upper tail where 
y = Lq + m, ip y = An [Qn] L ipi — > 0, thus it is 
also at the lower edge. 

2. If (k x , e) is a solution of N y + 1 = qL at the upper 

0, then by eqn. 



A20 



{k x 



edge, [Qn] {k x ) 

mcj), e) will be a solution of incommensurate Ny t 
1 = Lq + m. It is also at the upper edge because 



V'y 



in the lower tail. 



thus for N y — > 00, an increment of N y by 1 leaves the 
lower edge spectrum unchanged, while shifting the upper 
edge spectrum in k x by —2irp/q. 



Appendix B: Triangular lattice 

The Hofstadter model on a triangular lattice can be 
obtained by adding in each square plaquette a diagonal 
bond along the x — y direction with half-odd-integer vec- 
tor potential (in units of (f> = 2Ttp/q), such that the flux 
per triangle is <fi/2. Instead of eqn. [3j the Hamiltonian 




2 3 4 5 
k x (in 2w/q) 




0.5 




012345678 
k x (in 2-rr/q) 

FIG. 12: (Color online) Energy and entanglement spectrum 
on a triangular lattice with open boundary condition. Pa- 
rameters used are p/q = 5/8, N y = 32 and Ef = 1. In the 
energy spectrum (top panel), black dots represent bulk lev- 
els, red lines represent edge modes localized along the lower 
edge {y = 1), and blue lines represent those localized along 
the upper edge {y — N y — 32). Vertical gray line indicates 
the k x value at which Ef intersects the lower edge state. The 
entanglement occupancies (bottom panel) are computed for 
the lower half of the cylinder (1 < y < 16), color and symbol- 
coded according to a. The sudden color change happens when 
the lower edge mode crosses Ef. 



matrix is now 



H(k Xl N y ,z) = - 



( Ci 


Vl 





• • ■ ZV N y 




-l\ 


V*l 




V2 













1 


\ ZV N y 







V *Ny-l C Ny 


/ 



(Bl) 



where 



c y = 2 cos(k x + y(f>) , v y = 1 + 1' e -^ + ^ + 5<« , 

(B2) 
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£1(2/1,2/2), « = 1 



2/1 -i 



16 - 



24 



tlrj 
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-82(2/1,2/2), k = 1 
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16 
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FIG. 13: (Color online) Full system band projectors for p/q — 5/8, N y = 32, with periodic boundary conditions in y and 
k x — 2nn/q, Bj and Q v are projectors for the j th band and v lowest bands, repectively, see also Fig. [(j] Here we only plot the 
24 x 24 submatrix belonging to the first three magnetic unit cells. Since the Hamiltonian eqn. |Bl| can no longer be made purely 
real, the projectors are in general complex, so we only plot the real part of their matrix elements. Their signs are represented 
by color, red for positive and blue for negative, and their magnitudes represented by intensity. The individual bands B\ and 
B2 are better localized at k — 1 where the gap between them is maximal, while their sum, Q2, is better localized at k = 1.5 
where the gap between them is minimal. Note that at both k, the off-diagonal subblocks of both Bj within each 8x8 block 
tend to cancel (they have different colors). The adiabatic evolution of the wave packets are obvious: for example, at k — 1.5, 
the wave packets of Q2, in each unit cell, are at y = 3 and 6 (mod q. Same below). At k = 2, the wave packet at y = 3 is in 
the progress of moving toward y = 1 while the one at y = 6 is "frozen". At k — 2.5, the first wave packet arrives at y = 1. In 
the next « — > n + 1 sub-cycle, the first wave packet will be frozen and the second one will migrate in the diagonal line by —2 
(negative of the total Chern number) . 



with t' = 1 for triangular lattice, and for square lattice. 
The discussion of edge spectrum in Appendix |A| remains 
essentially the same, except the transfer matrix, eqn. |A3| 
now becomes 



M„ 



V 



1 



/ 



(B3) 



In the pathological case where certain Vy = 0, the open 
edge Hamiltonian (i.e., z = 0) reduces to two blocks, 



{1, . . . , y} and {y + 1, . . . , N y }, each of which can be in- 
dividually solved; alternatively one can shift t slightly 
away from 1. Note that while detM y is no longer one, 
the (/-step transfer matrix Q still has unimodular deter- 
minant, 



|det Q\ 



v o v i 



■ v 



q-l 



viv 2 ■■■v q 



1 



where we used vq — v q . Consequently 

IQ22I = IQ11I" 1 , 



(B4) 



(B5) 
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and eqns. |Al8] and pO0| hold up to a phase The con- 



clusion thus remains unchanged that lower edge states 
are unchanged while upper edge states shift in k x with 
incommensurate N y . 

In Fig. |12[ we plot the open boundary energy spec- 
trum, and its entanglement spectrum with v = 5 filled 
bulk bands for p/q = 5/8 on the triangular lattice. The 
Chern numbers of individual bands are either C = —3 
or 5, which are equivalent modulo q = 8, and the lowest 
band has C — —3. Both are in agreement with what 
we observed in the square lattice case, namely, all band 
Chern numbers are equivalent modulo q, and that the 
lowest band Chern number is the one with smaller mag- 
nitude. The number of edge spectral flows in each gap 
is the total band Chern number below the gap, and the 
entanglement spectral flow mimics the behavior of edge 
spectral flow, and has an index discontinuity at k x where 
Fermi energy intersects the lower edge mode. Note that 



for v = 5 filled bands, the total Chern number is 1, as re- 
flected in the number of edge and entanglement spectral 
flows. This agrees with our observation that the total 
Chern number of p filled bands is one, see discussion fol- 



lowing eqn. 20 in the text 



As in the square lattice case, the band projectors and 
their sums also flow under adiabatic k x pumping, and 
the number of wave packets crossing any given bound- 
ary during one cycle of the pumping reflects the Chern 
number of the projectors. Fig. [13] shows the flow of the 
lowest two band projectors, B\ and B2, and their sum Q2, 
at k = K- 2n/q with k = 1 and 1.5. Both Bj and Q v have 
better localization at either integer or half-odd-integer k 
where its gap from neighboring bands are maximal. 

We thus conclude that the observations as detailed in 
the text using square lattice are robust and insensitive to 
the underlying lattice used. 
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